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We consider the nonlinear Schrodinger equation with a point-like source term. 
The soliton interaction with such a singular potential yields a critical solution 
behavior. That is, for the given value of the potential strength and the soliton 
amplitude, there exists a critical velocity of the initial soliton solution, around 
which the solution is either trapped by or transmitted through the potential. In 
this paper, we propose an efficient method for finding such a critical velocity by 
using the generalized polynomial chaos method. For the proposed method, we 
assume that the soliton velocity is a random variable and expand the solution 
in the random space using the orthogonal polynomials. The proposed method 
finds the critical velocity accurately with spectral convergence. Thus the com- 
putational complexity is much reduced. Numerical results for the smaller and 
higher values of the potential strength confirm the spectral convergence of the 
proposed method. 

Keywords: Nonlinear Schrodinger equation. Singular potential. Generalized 
polynomial chaos. Stochastic collocation method. Split step Fourier method, 
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1. Introduction 

Nonlinear Schrodinger equation (NLSE) describes a broad range of physical 
phenomena, e.g. nonlinear modulation of collisionless plasma waves [1], self 



trapping of a light beam in a color dispersive system [2'] , helical motion in a very 
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thin vortex filament , propagation of heat pulses in an-harmonic crystals , 
modulation instability in water waves [lo| . etc. In optical fibers, the soliton 
solutions of the NLSE provide a secure means to carry bits of information over 
many thousands of miles ■ Termed as the Gross-Pitaveskii equation, the NLSE 
with an appropriate potential can be utilized to describe the dynamics of the 
Bose-Einstein condensate, both with the attractive and repulsive nonlinearities 
is our objective in this paper to solve the Gross-Pitaveskii equation 



equipped with a point-like potential to find the critical values of the soliton 
velocities when the amplitude of the point-like potential is either very small 
{r^ 10^^ — 10^^) or large (~ 2.5 — 4.5) compared to the soliton amplitude which 
is the unity in our paper. 

We know the soliton solution of the homogeneous NLSE 

idtu+^dlu + u\u\'^ ^0, -oo < X < oo, t > 0, (1) 

with initial condition uq given by 

uoix) = Asech {A (x)) e(*'^+*^=^) (2) 

is given by 



u{x, t) — Asech {A {x — Vt)) exp I icj) + 



A > 0, V eR, (3) 



where A is the soliton amplitude, V the soliton velocity, and (f> the phase lag. 
Consider a perturbed NLSE, that is, the Gross-Pitaveskii equation by adding 
an external potential, —eS{x)u, 

idtu + ^d'^u + ulul"^ = —ed{x)u, 
u{x,0) = uo(x), 

where 6{x) is the Dirac delta function with a constant e G R. Such an ex- 
ternal potential represents the impurity or defect in the optical fiber. The 
well-posedness of the equation 

idfU + ^d^u + u\u\P^^ = —e6{x)u, 
u{x,0) = uo{x), 

with p > 1 and initial data uq in _ff^(]R), has been extensively studied and is 
based on the knowledge of the self-adjoint (in L^) operator —dxx + e<5. Using 
[3] , Le Coz et al. proved the existence of a time T > and of a unique solution 
to Eq. H (where e G M) in C{[0,T), H\R)) n ([0, T), iJ-i(R)) satisfying 
limf^T- ||9j;u||2 = oo. Moreover the energy is conserved in time. This result 
was extended to p > 1 by Fukuizumi et al. in Q. For p — 3 (more generally 
p £ (1, 5)), global existence in also holds by Gagliardo-Nirenberg's inequality 
and energy conservation. Global existence in H^, is also discussed by Goodman 
et al. using a fixed-point argument and time-invariance of the L^-norm and 
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of the Hamiltonian derived from the NLSE. Notice that the study of stabihty 
of nonhnear bound states which are solutions of the form e'xjp(—±ujt)(f)i^{x) with 
w > 0, and for which: 

- ■^dxx(j)u, - e(j)uj - \(l)u?(f)i^ = ^^4>uJ (6) 

plays an important role in the theory of NLSE with defect and could possibly 
be useful numerically. Explicit formulas and stability analysis for 0^ can be 
found in 



If we now take a soliton approaching the impurity from the left as an initial 
condition uq: 

uo{x) = Asech {A {x - xo)) e(''^+'^^\ xo < 0, (7) 

then until the time io = f?, the solution will still be given by Eq. |31 In this 
paper we consider A = \ and (j) = Q. Thus the soliton velocity V and the 
strength of the impurity e are the only parameters of the problem. 

For io > f?'; the effects of the potential are highly visible and a lot of 
research has been done on the transmission and reflection coefficients of the 
^-potential by the standard scattering theory [l2j . Malomed and his co-workers 
(si [l7j showed mainly numerically, that for any given velocity V {> 0), there 
exists a threshold value ethr {> 0) of e, for which the soliton can marginally pass 
through the defect. So for the given velocity if e < Cthn the soliton can pass 
through the defect and the soliton gets trapped otherwise. They considered 
the soliton-soliton collisions within the coupled NLSE. In the limiting condition 
one soliton has very large amplitude and is very narrow accordingly, while the 
soliton governed by the other equation has finite amplitude and width. In this 
limiting condition the two coupled NLSE are reduced to a single equation, in 
which the narrow soliton will be represented by the (5-function, 

idtu+]^dlu + u\u\'^ = -e5{x)u. (8) 

Holmer and his co-workers studied the NLSE with V ^ \ and V <C 
1, e ^ 1 ^lij. They showed for high there exits the bound state which is 
given by 

u{x,t) = e'-^^sAsech {\\x\ + tanh"^ (e/A)) , < A < e, 

and this bound state is "left behind" after the interaction (see bottom right 
figure of Figure [3]). Also they proved in [l^ that ioi V ^ 1 and e <C 1, the 
solution can be approximated by the soliton solution of the homogeneous NLSE 
(e = 0). To solve Eq. [S] for any given e (> 0) and V {> 0), we consider three 
cases: (a) small value of e, where e < 0.3 (b) moderate value of e, where 0.3 < 
e < 3.5 [3| and (c) large value of e, where e > 3.5. For solving Eq. [51 one can 
use the Split Step Fourier Method (SSFM) to reduce the computational time. 
To get tthr for any given V with certain accuracy one must conduct a series 
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of simulations. The number of simulations increase with the increase of the 
level of accuracy. In addition to conduct a series of simulations with small time 
steps, one needs a large amount of the computational time. This is our main 
motivation to propose a suitable method to overcome such a high computational 
complexity by using the generalized polynomial chaos (gPC) metho ds Jl | 



The gPC method belongs to the class of non-sampling methods [20|, l21| . In 
this method the stochastic quantities are expanded by orthogonal polynomials. 
Different types of orthogonal polynomials can be chosen for better convergence. 
The gPC expansion is a spectral representation in random space and exhibits 
fast convergence when the expanded function depends smoothly on the random 



parameters When the gPC method is applied to solve any differential 



equation, the main computational work is needed to solve the expansion co- 
efficients of the gPC expansion. A common approach is the Galerkin method 
that minimizes the residue in the polynomial space. The stochastic Galerkin 
(SO) approach, however, would be extremely difficult to use when the governing 
stochastic equations take complicated forms. In our case, the NLSE contains 
the nonlinear term |upit. For the SG method, it is very hard to get the corre- 
sponding explicit deterministic equations after expanding the nonlinear terms. 
So that, in this work we use the high-order stochastic collocation (SC) approach 
(iot that combines the advantages of both the Monte Carlo sampling and the 
gPC-Galerkin methods. The gPC method reduces the number of simulations 
for finding the critical velocity, Vc, for any given value of e thanks to the high- 
order convergence of the method. Since the equation has only two parameters, 
i.e. e and V, we treat at least one of them as a stochastic variable in the gPC 
framework. In the present work we consider V as the stochastic variable and 
let e be fixed. So for any given e, we find Vc, the critical value of V around 
which the soliton is either transmitted or trapped. Thus it is obvious that for 
V > Vc, the soliton passes through the defect. By adopting this idea we develop 
a step-by-step gPC collocation method to find the critical velocity of the soliton. 

In I3| the relation between ethr and V was obtained only for the moderate 
values of e, i.e. for those comparable to the soliton amplitude A — 1. But the 
results of the numerical simulations for very small or large values of V were not 
obtained, perhaps due to the huge computational burden. By the gPC method, 
we were able to reduce the overhead computational time, for having detailed 
simulations performed for large and small values of e to find the corresponding 
critical velocity Vc- 

Since the analysis for the moderate values of e are already done Q, we do 
not intended to repeat the analysis for those values of e in this paper. Here we 
mainly focus on the small and high values of e. For the small values of e, the 
gPC takes much longer time than the gPC method for the large values of e due 
to the extremely small critical velocities. 

This paper is organized as follows. In Section 2, we discuss the SSFM. 
Section 3 describes the gPC collocation method. Section 4 contains the gPC 
collocation algorithm for the NLSE with the singular potential term to detect 
the critical velocity for the given value of e. Section 5 presents the numerical 
results. Concluding remarks and future works are presented in Section 6. 
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2. Split Step Fourier Method 

The SSFM is a pseudo-spectral numerical method used to solve nonlinear 
PDEs like the NLSE. Eq. [T]can be rewritten as 

^ = ^[N + D]u, (9) 

where D — \-§^ and N — The solution of Eq. [3] can be written as 

where u{x^ 0) is the initial condition. Since D and N are the operators, they do 
not necessarily commute. However the Baker-Hausdorff formula can be applied 
to show that the error will be of order dt^ if we are taking a small but finite 
time step dt We therefore can write 

u(x,t + dt)Ke"^'^e"^'^u{x,t). (10) 

The part of this equation involving N can be computed directly using the wave 
function u{x, t) at time t. To compute the exponential involving D we use the 
fact that in the frequency domain, the partial derivative operator ^ is converted 
into «fc, where k is the frequency associated with the Fourier transform. Then 
we take the Fourier transform of u{x, t) recover the associate wave number, and 
compute 

e-¥<n^^W[u{x,t)] , 

where F denotes the Fourier transform. Then we take the inverse Fourier trans- 
form of the expression to find the solution in the physical space, yielding the 
final expression 

We apply SSFM to Eq. |8] where the nonlinear operator N — |up and the linear 
operator L = \-§^ + €5[x). In our numerical simulations we use the high-order 
SSFM, such as the Strang splitting based on: 

gZ(L+JV)At ^ ^,L^^,NAt^^L^ ^ 0{M''{\L, [L, N]] + [N , [TV, L]])) . 

where [L, N] — LN — NL denotes the commutator between L and N. Thus, 
from t to t + At 

u{x,t + At) ^ e^'^^+^^^'uix.t), 

« e'^^e'^^*e'^^u{x,t). (11) 
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3. gPC collocation method 



We solve Eq. 0] with the initial condition given by Eq. [7] for both small and 
large values of e by the gPC collocation method. We use the gPC method for 



the solution of the NLSE using the Wiener- Askey scheme [19|, |2l[ , in which Her- 
mite, Legendre, Laguerre, Jacobi and generalized Laguerre orthogonal polyno- 
mials are used for modeling the effect of continuous random variables described 
by the normal, uniform, exponential, beta and gamma probability distribution 
functions (PDFs) , respectively [20| . These orthogonal polynomials are op- 
timal for those PDFs since the weight function in the inner product and its 
support range correspond to the PDFs for those continuous distributions. 

Following the standard gPC expansion, we assume that u{x,t,^) is suffi- 
ciently smooth in ^ and has a converging expansion of the form 

oo 

u{x,t,^) ^^Uk{x,t)Pk{£,), 

fe=0 

where the orthonormal polynomials Pk (^) correspond to the PDF of the random 
variable ^ and satisfy the following orthogonality relation: 



E[PkPi] := J PkiOPiiOpiOd^^Ski. 



Here 6ki is the Kronecker delta and p{£,) is the weight function. Note that the 
polynomials are normalized. 

For the stochastic collocational approach we approximate Uk{x,t) as, 

Q 

Uk{x,t)^^u{x,t,p>)Pk{p')a,, fc = 0,---,Q, (12) 

where Q -I- 1 is the total number of the collocation nodes. Here {p' , a-' } is a set 
of nodes and weights, where p' and denote the j-th node and its associated 
weights, respectively, in the random space F such that 



is an approximation of the integral 

I[f]^ J^fip)pip)dp^E[fip)], (14) 

for sufficiently smooth functions f{p), i.e, 

[/] ^ / [/] , Q ^ oo. 

In this paper we consider V as the stochastic variable and we choose a 
collocation nodal set { V^-' , a-' } in space F, where are the jth collocation 
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Figure 1: Soliton interaction with the defect with small strength (e = 0.1), where initial 
velocity of the soliton is zero. The soliton is trapped and an oscillatory movement is observed. 

points and a-' the corresponding weights. For each j — 0, ■ ■ ■ ,Q, we solve the 
problem given by Eqs. H] and [7] with the parameters e and and let the 
solution set be {uq, ■ ■ ■ ,uq} where uj is the solution for V = Vj. For solving 
this deterministic equation, we employ the high-order SSFM. The approximate 
gPC expansion coefficients are 

Q 

Um{x,t) = {x,t,Vj) (pm (^jO^j, TO = 0, • • • , Q, 

where {(pm) are the orthonormal polynomials. And finally we construct the Qth 
order gPC approximation 

Q 

u{x,t;V) w ^ Um{x,t)(j)m {V) , whcrc ^ = {Vq,V2, ■ ■ ■ ,Vq} . 

4. gPC collocation algorithm for solving NLSE 

The following algorithm describes how to calculate the critical velocity by 
using the gPC collocation method. 

We use the gPC method to find the critical velocity Vc efficiently for any 
given e. Here the soliton velocity V is the stochastic variable. Suppose we know 
in advance that the critical velocity Vc hes between Va and Vb {Va < Vt) and 
consider V has a uniform distribution over [Va, Vf,]. Since the distribution is 
uniform, we use the Legendre polynomials for expanding the solution in the 
random space. For this purpose we choose + 1 Gauss-Legendre quadrature 
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Figure 2: Top: Soliton solution witiiout any defect. Middle: Soliton passes through the defect 
with the initial velocity V = 0.001 and the defect amplitude e = 0.1. Bottom: Soliton trapped 
by the defect with the initial velocity V = 0.003 and the defect amplitude e = 0.5 
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Figure 3: Top left: Soliton transmitted tlirougli the defect when e = 0.08 and V = 5 X 10~^. 
Top right: Sohton transmitted through the defect when e = 0.1 and V = 8 X 10~^. Bottom 
left: Soliton is trapped by the defect when e = 4.5 and V = 0.220048. Bottom right: Soliton 
is transmitted through the defect when e = 4.5 and V = 0.23995187. For e = 4.5, radiation 
effect is clearly visible. 



points with the weights. Let the set {a^, i^iji^o describe the {N + 1) quadrature 
points ai and the corresponding weights w^. 

Now find the solution of Eq. |4]for each V = aihy using the high-order SSFM. 
For this purpose one must use a sufficiently large computational domain and 
sufficiently long time interval. We set up the domain size and the computational 
time in such a way that no solution leaves the domain yet with the given final 
time. For example when e = 0.3, we use the domain size [—L, L] — [—40, 40] 
and the final time tf = 12000. We are solving the NLSE for Uj {x, t, Vj) for all 
Vj with the same final time. 



We reconstruct the soliton solution for each simulation for x e 



'L, L 



at 



the final time. L is chosen in such a way that only the trapped solutions exist 



inside 



-i, L 



. We know if the solution is trapped, it would stay around the 



position of the defect (in our case at a; = 0). So L must be close to zero. In 
our computation, we choose L where the mean solution vanishes near x — 0'^ . 
For the i-th quadrature point ai, we denote the solution by Ui {x,Tf,ai). 
Evaluate the approximate gPC expansion coefficients by 



Q 

Urn {x,Tf) =^Wi {x,Tf,ai)Lm {cti) UJi, 




Figure 4: The nonlinear interaction of the soliton with the defect (the dotted hne). Top: 
The nonlinear interaction is prominent when the soliton hits the deffect (e = 0.3) with small 
velocity (~ 10~^) compared to the interaction with the high velocity (~ 10~^) where € = 4.5 
(middle). Bottom: The interaction of slowly moving soliton (V ~ 10"^)) with the defect with 
high value of e (= 4.5). 
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where {im}„j=o Legendre polynomials and uJi are the quadrature 

weights. The full gPC solution is given by 

Q 

u{x,Tf,a) ^'^Uk{x,Tf) Lk{a). (15) 

fc=0 

The mean solution is given by the 1st mode jl^l: i-e 

Q 

uq {x, Tf) = ^ Ml {x, Tf, ai) Lq (a;) oJi. (16) 

i=0 

From Eq. 1161 one can construct the average energy E of the system between 
—L, L and [—L, L] at the final time, that is, 

El^-J \uo{x,Tf)\^ dx, E^, = - j \uQ{x,Tf)\^ dx. (17) 

Suppose that among N solutions, A^i solutions are trapped inside — i, L 
Then A^i can be estimated for large — !> oo by 

Ni^E^ Ni^ Vc-Vg 
N El' N Vb^Va' 

where Vc is the critical velocity for given e. So Vc is evaluated by 

V, = Va + {Vb-Va)^. (18) 

If we increase the number of quadrature points, then the critical velocity can be 
determined more accurately. For our simulations we used 24 Gauss Legendre 
quadrature points and obtained spectral accuracy of ^ 10~^^. Figure [7] shows 
the spectral convergence of the error of the critical velocities with the increasing 
number of the quadrature points. 
Remark: 

The solution u{x, V) has possibly a jump atV^Vc for t ^ oo because of the 
critical behavior of the soliton solution around the potential. This means that 
the spectral reconstruction of u{x,t,V) for any V G [Va,Vi,] using ui{x,t),l = 
0, • • • , Q may fail to converge to the right solution due to the discontinuity at 
V = Vc- This was also addressed in our previous work for the critical behavior of 
the soliton solution for the sine- Gordon equation 0/. Here note that the proposed 
method in this paper uses only the first moment uo{x,t) to estimate the critical 
velocity Vc but not the reconstruction of u(x,t,V) . The mean solution, uo{x,t) 
is convergent. 

In Eq. [THl the convergence of Vc mainly depends on i? As the defi- 

nition in Eq. 1171 the convergence of R then depends on how UQ{x,t) converges 
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Figure 5: First mode (mean) of the gPC expansion for different e. Top and middle: The 
Legendre chaos. Bottom: The Hermite chaos. Right figures of the top and middle panels 
show the locations of L for different e. 



with N. In our previous work 15[, it was proven that uq{x) converges fast 
enough although the original function u{x, V) is discontinuous in the random 
variable V . As we will discuss in the next section, numerical results in Section 
5 (Figure 7) implies that R shows spectral convergence with N . 



5. Numerical results 

We first consider the high value of e, say e — 2.7. By doing few Montc- 
Carlo simulations we roughly estimate the interval F G [K, VJ,] , K G [K, Vb] 
where Vc, the critical velocity may be located. For e — 2.7, we use Va = 0.1 
and Vb = 0.14. Since for moderate and high values of e, the simulation time is 
relatively less than the simulation time with smaller range of e, we follow the 
same procedure to find the suitable intervals. But for the small value of e, i.e. 



12 




Figure 6: Critical Velocity vs. e where e € [0.05, 4.5]. 



e < 1.0, where the simulation time is long, we use the extrapolation of Vc from 
the previous e to get the rough estimate of the interval. 

To apply the gPC collocation method, one also needs to find the value of L . 
We do not have any fixed L which can serve for all e. Instead, we have different 
L for different e. A heuristic approach is used to find L . For the given value of 
e, we construct the mean solution by Eq. \Wi Since some solutions are trapped 
and some of them are transmitted, there are few bumps near the defect and few 
bumps are far from the defect. Clearly there exists a separation point between 
these two groups of bumps. Ideally the x— coordinate of this point would be 
zero but due to the domain truncation, radiation effect etc. it may not be equal 
to zero. By observing the graph carefully we can easily find the separation 
point which we use as L . For the small and moderate values of e, determining 
accurately L is easy, but for the high values of e, we need extra care. For 
the high value of e, the values of Va, Vb are also high and we can not run the 
simulations for a long time because some solutions may leave the domain and 
re-enter the domain from the other side due to the periodic boundary conditions. 
So in this case we need to study the bumps carefully to locate L . In Figure [5l 
the zoomed graphs of the mean solution of each e are given in the right panel 
of the top and middle figures. We find that for e = 0.3, L = 12 and for 
g = 0.5, L = 13.5. Similarly for e = 2.7, L = 10 and for e = 3.0, L = 15. 

Figure [T] presents the interaction of the soliton with the (5-function. Here 
we choose the initial velocity, Vq = and the potential strength e = 0.1. The 
soliton is located at xq — —0.3 initially, which is inside the influence zone of 
the potential. The nonlinear interaction is observed and the soliton solution 
exhibits an oscillatory behavior along the line a; = 0. This case was discussed 
in 13, 3|- But such an initial condition may not necessarily satisfy the given 
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equation. The initial position of the soliton must be out of the influence zone of 
the potential and the soliton must be allowed to move freely before it hits the 
defect. In all cases we consider the starting point of the soliton (xq) is far from 
the position of 5-function, i.e. outside the influence region of the potential. 
Figure [5] shows the behavior of the soliton solutions in three different cases. 
When e = 0, that is the case when there is no (5- function, the soliton solution 
passes unperturbedly. But for nonzero e, the soliton behaviour depends on its 
initial velocity. For e = 0.1, the soliton passes through the defect for V — 0.001 
and for e = 0.5 and V — 0.003, soliton is trapped by the defect. For both cases, 
the soliton passed or trapped as a whole. There is no radiation due to the small 
soliton velocities 

Figure [3] represents the long time simulations for (e, V) = (0.08, 5 x 10^^) 
(top left), (0.1, 8 X 10-5) (top right), (4.5, 0.220048) (bottom left) and (4.5, 0.23995187) 
(bottom right). For the case that e is small and V is also very small accordingly, 
the soliton is transmitted through the defect without any radiation. But for the 
high value of e, usually greater than 2.7, where the critical velocity is also high, 
the radiation effect is observed due to the soliton-defect interaction. The bot- 
tom panel of Figure [3] exhibits the radiation effect for e = 4.5. For both the 
"trapped" and "transmitted" situations, the radiation effect is observed. The 
bound state effect is also observed in the bottom right, the details of which was 
discussed in [l^ . 

Figured] shows the nonlinear interactions of the soliton with different soliton 
velocities. When the soliton velocity is small, nonlinear property dominates as 
shown in Figure [3l During the time of interaction with the defect (the dotted 
line), the soliton velocity increases and after crossing the defect, the velocity 
turns into its previous value. When the soliton velocity is high, the linear effect 
dominates and the soliton velocity does not changes during the collision but the 
direction of the propagation changes. That is the soliton continues its motion 
with the same velocity. When a slowly moving soliton hits the defect with high 
strength (e = 4.5), the soliton is trapped by the defect but due to the nonlinear 
interactions, radiations and transmissions are also seen (bottom figure). 

Figure [5] shows the mean solutions at the final time. This is the first mode 
of the solution by the gPC collocation method. Here we used F as a stochastic 
variable, V G [Va, Vb] and T4 and Vb are different for different values of e. We 
used both the Legendre and Hermite chaos. We need to consider the uniform 
distribution and normal distribution for the Legendre and Hermite chaos respec- 
tively. In Figure [HI the figures in the top panel are obtained using the Legendre 
chaos for e = 0.3, 0.5. Those solitons that are trapped by the defect are con- 
fined around the position of the defect. In our case, the defect, the (5-function 
is located at a; = 0. There are multiple peaks in the mean solution, but around 
a; = the peaks are higher than the others, which implies that some solitons 
are trapped, and the rest are transmitted. These figures are used to locate the 
position of L . If we see the zoomed figure in the right panel, we easily locate 
L for different e. 

For the middle panel figures in Figure [5l we plotted the mean solutions and 
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zoomed one for e = 2.7, 3.0 and 4.5. The sharp peaks at x = imply that 
the most of the solutions are trapped in that range of V and some of them are 
transmitted. We already mentioned that in this region of such a large value 
of e, the radiation effects are visible, which are also showed in the figure. The 
values of L are pointed for different e values in the figure. Same explanation 
for e = 0.5. 

Next we consider the case that V is normally distributed and we use Her- 
mite polynomials [l^ and the Gauss- Hermite quadrature points 11 1 . Let Va = 
a, Vb — f3 and V G [a, /?], ^ G [—1, 1], 7 G (—00, cx)). The linear transforma- 
tion between V and ^ is given by 



and the transformation between ^ and 7 is given by 

7 = e^o, 



= 0, 



Or we have, 



27 



Thus we have, 



= 0, 



I3~a 



7 = 0. 



-1 + ^1 + 472 
27 



+ -(« + /?), 



(19) 



where 7 has the normal distribution with mean and the standard deviation 
(SD) 0.1. For the simulation we consider e = 0.3 and y - iV [0, 0.1]. The 
figure in the bottom panel of Figure [5] shows the mean solution at the final time 
obtained by the Hermite chaos. Although the mean solutions obtained from the 
Legendre and Hermite chaos arc different, we observe that the location of L is 
same for both cases. 

Using a series of those simulations above for different values of e where 
e e [0.05, 4.5], we determine the critical velocities with respect to different e. 
The results are plotted in semi- logarithmic scale in Figure |6l It is observed that 
for the small values of e where e < 0.1, the curve is very stiff and the slope 
changes sharply around e — 0.1. From e > 0.1, the curve increases steadily. 
The "trapped" and the "untrapped" regions are clearly shown in the figure. The 
V — e graph is the boundary of those two regions. 

5.1. Convergence analysis 

We define the error of the critical velocities by 



Error'(7V) = \V^iN) - V^'{N - 1)| , 
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Table 1: Convergence of Vc with for the Legendre Chaos, e = 0.3, 1.0 and 4.5. 



N 




Vc X 10'' 






e = 0.3 


e = 1.0 


e = 4.5 


2 


1.658675134594813 


23.83399810435849 


233.3998104358486 


4 


1.786459011090171 


24.00997282851472 


236.0392987896415 


8 


1.788091882999389 


24.01306953886891 


236.7359886666223 


12 


1.788112748469627 


24.01312403314395 


236.9198791482719 


16 


1.788113015096694 


24.01312499210545 


236.9684168267407 


20 


1.788113018503758 


24.01312500898075 


236.9812282904446 


24 


1.788113018547295 


24.01312500927771 


236.9846098613687 



where N is the number of collocation points. Figure [7] shows the convergence 
of errors obtained by the Legendre and Hermite chaos. We do the convergence 
analysis for various values of e. We choose e — 0.3 (small) , e = 1.0 (moderate) 
and e = 4.5 (high). For the Legendre chaos, the critical velocities for different TV 
are presented in Table 1. For e = 0.3 and e = 1.0, we calculate the errors for both 
the Legendre and Hermite chaos and for e — 4.5 we use the Legendre chaos. For 
Hermite chaos, we expect to have the similar results. The graphs are plotted in 
semi-logarithmic scale. Figure [7] shows all the graphs are a straight line, which 
confirms spectral convergence but the convergence rates are different for different 
cases. For e = 0.3 and e = 1.0, Hermite chaos exhibits slower convergence rate 
than the Legendre chaos. Also if we compare the graphs for the Legendre 
chaos for different cases, it is found that the convergence rate decreases with 
the increases of the value of e. That is, the smaller is the value of e, the faster 
convergence is obtained. One of the possible reasons is because of the radiation 
effect. As e increases, the radiation effect becomes visible and it makes difficult 
to locate the position of L accurately. According to our numerical results, our 
main result is stated by the following; The numerical scheme stated in Section 
4 to find the critical velocity Vc has the spectral convergence and the rate of 
convergence decreases with increase of the value of e. 

6. Conclusion 

In this paper we studied the NLSE with the singular potential. We proposed 
an efficient method of determining the critical soliton velocities, Vc, by using 
the gPC collocation method. We studied the wide range of e, i.e. e G [0.05, 4.5]. 
For e < 0.05 the numerical simulations demand a huge computational time due 
to the very small soliton velocity (V ^ 10~^°). We studied the convergence 
analysis to prove the merit of our proposed numerical scheme. We found the 
spectral convergence in all cases. The main development of this paper is the use 
of the gPC collocation method to determine the critical velocity of the soliton 
for given e with the desired level of accuracy. We obtained Vc accurately with 
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Figure 7: Spectral convergence of the critical velocities for e = 0.3,0.1 and 4.5. Graph shows 
the spectral convergence for both the Legendre and Hermite chaos. Note that the Legendre 
chaos shows faster convergence than Hermite chaos. 



a small number of simulations. In our future work, we will further study the 
case that e ^ 0.05. Also for the high values of e, where radiation effect is 
prominent and the convergence of the proposed method becomes slower due to 
the radiation effect, an efficient numerical method dealing with this effect will 
be investigated. 

Acknowledgement: The first author is grateful to Gino Biondini for de- 
veloping and implementing high-order SSFM. 
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